Baseline dynamics

\[\begin{align} S' & = \mu N^\star - (\lambda + \mu) S \\ E' & = \lambda S - (\kappa E + \mu E) \\ I_S' & = p \kappa E - (\gamma_S + \delta_H + \mu_{I_S} + \mu) I_S \\ I_A' &= (1 - p) \kappa E - (\gamma_A + \mu) I_A \\ H' &= \delta_H I_S - (\gamma_H + \mu_H + \mu) H \\ R' & = \gamma_S I_S + \gamma_A I_A + \gamma_H H - \mu R \\ D' &= \mu_{I_S} I_S + \mu_H H \\ \lambda &:= \frac{\beta_A I_A + \beta_S I_S}{N^{\star}} \\ N^{\star}(t) &= S + E + I_S + I_A + H + R \end{align}\]

Stan ODE notation beginning model

state Stan Initial Condition Stan Postulated
\(S\) y[1] \(S_0\) y_init[1] \(N - (I_{S_0} + I_{A_0} + E_0)\)
\(E\) y[2] \(E_0\) y_init[2] Unif(74, 2100)
\(I_S\) y[3] \(I_{s_0}\) y_init[3] 22 (cdmx data)
\(I_A\) y[4] \(I_{a_0}\) y_init[4] Unif(74, 2100)
\(H\) y[5] \(H_0\) y_init[5] 0
\(R\) y[6] \(R_0\) y_init[6] 0
\(D\) y[7] \(D_0\) y_init[7] 0
\(I_s^{[cumulative]}\) y[8] \(I_{s_0}^{[cumulative]}\) y_init[8] 74 (cdmx data)

Observation Model

\[\begin{align} Y_t & \sim \text{Poisson}(\lambda_t), \qquad \\ \lambda_t & = \int_0^t \rho \delta_e E \\ & p \sim \text{Uniform(0.3 0.8)} \\ & \kappa: \sim \text{Gamma(10,50)} \end{align}\]

March 10 – March 30 parameters estimation

stan notation symbol Prior Fixed
theta[1] \(\beta_S\) Normal(1, 0.3)
theta[2] \(\beta_A\) Normal(1, 0.3)
theta[3] \(\kappa\) Gamma(10, 100)
theta[4] \(p\) Unif(0.3, 0.8)
theta[5] \(\delta_H\) Gamma(10, 40)
theta[6] \(\mu_{I_S}\) Gamma(10, 80)
theta[7] \(\mu_H\) Gamma(10, 350)
theta[8] \(\gamma_S\) Gamma(10, 350)
theta[9] \(\gamma_A\) Gamma(10, 350)
theta[10] \(\gamma_H\) Gamma(10, 350)
y_init[2] \(E_0\) Unif(74, 2100)
y_init[3] \(I_S\) ———– 22
y_init[4] \(I_A\) Unif(74, 2100)

Postulated densities

## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
mu_ <- function(lambda, ylabel= "mu_h"){
    # We suppose mu_. as a exponential r.v.
    x <- seq(from = 0, to = 1, by=.001)
    mu_h_ <- dexp(x, rate = lambda)
    df <- data.frame(x=x, mu_h = mu_h_)
    mean_ <- 1/lambda
    mean_time <- 1/mean_ 
    text <- paste(mean_time,"days")
        p <- ggplot(data=df,
                aes(x = x,
                    y = mu_h)) + 
        xlab("x")+
        ylab(ylabel) +
        geom_line(colour = "#000000") +
        geom_area(alpha = 0.6, fill = "lightgray") + 
        geom_vline(xintercept = mean_, linetype="dotted", 
                color = "red", size=0.5)+
        annotate(geom = "text", x = 1.30 * mean_, y=0, label = text,
              color="black")
    return(p)
}
#
gamma_ <- function(alpha, beta, ylabel="gamma_S"){
    # We suppose mu_. as a exponential r.v.
    theta <- 1 / beta 
    x <- seq(from = 0, to = 1, by=.001)
    gamma_s_ <- dgamma(x, shape = alpha, 
                        scale= theta)
    df <- data.frame(x=x, gamma_s = gamma_s_)
    mean_ <- alpha * theta
    mean_time <- 1/mean_ 
    text <- paste(mean_time,"days")
    p <- ggplot(data=df,
                aes(x = x,
                    y = gamma_s)) + 
        geom_line(colour = "#000000") +
        geom_area(alpha = 0.6, fill = "lightgray") +
        geom_vline(xintercept = mean_, linetype="dotted", 
                color = "red", size=0.5) + 
        annotate(geom="text", x=1.30 * mean_, y=0, label = text,
              color="black")
        return(p)
}
pi_mu_h <- mu_(lambda = 25, ylabel= "mu_h")
pi_mu_i_s <- mu_(lambda = 35, ylabel= "mu_i_s")
pi_gamma_s <- gamma_(alpha = 10, beta = 100,  ylabel= "gamma_S")
pi_gamma_a <- gamma_(alpha = 10, beta = 50,  ylabel= "gamma_A")
pi_gamma_h <- gamma_(alpha = 10, beta = 250,  ylabel= "gamma_H")

fig6 <- ggplotly(pi_mu_h)
fig7 <- ggplotly(pi_mu_i_s)
fig8 <- ggplotly(pi_gamma_s)
fig9 <- ggplotly(pi_gamma_a)
fig10 <- ggplotly(pi_gamma_h)
fig6

Stan estimation

theta mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta_s 8.690483e-01 0.0051487 0.0978857 6.821660e-01 8.046863e-01 8.683529e-01 9.370239e-01 1.061825e+00 361.44444 1.0014380
beta_a 7.738431e-01 0.0046027 0.1046186 5.690834e-01 7.016002e-01 7.722679e-01 8.464557e-01 9.786888e-01 516.63452 0.9995944
kappa 2.143749e-01 0.0041385 0.0453536 1.506806e-01 1.847692e-01 2.042473e-01 2.344915e-01 3.419625e-01 120.09807 1.0212692
p 4.387130e-01 0.0024290 0.0526256 3.089264e-01 4.087240e-01 4.528404e-01 4.803055e-01 4.989015e-01 469.38544 1.0011927
delta_H 4.603172e-01 0.0123912 0.1174493 1.971069e-01 3.909769e-01 4.616815e-01 5.355033e-01 6.810346e-01 89.84121 1.0140202
mu_{I_S} 2.602712e-01 0.2386425 1.3275922 4.162000e-03 2.803170e-02 5.691990e-02 1.042763e-01 3.605100e-01 30.94808 1.0442815
mu_H 2.327870e-01 0.2376622 1.3076330 2.533600e-03 1.305290e-02 3.280390e-02 6.643520e-02 3.525362e-01 30.27273 1.0449172
gamma_S 1.227483e-01 0.0034395 0.0412130 5.379330e-02 9.309530e-02 1.197828e-01 1.465723e-01 2.144921e-01 143.57880 1.0148258
gamma_A 5.106448e-01 0.0079512 0.0939779 3.163330e-01 4.584978e-01 5.113880e-01 5.662136e-01 6.887404e-01 139.69809 1.0283537
gamma_H 5.079869e-01 0.0089289 0.1639180 2.524056e-01 3.875971e-01 4.880279e-01 6.033939e-01 8.925461e-01 337.02411 1.0089564
R_zero 1.916769e+00 0.0988256 0.6028527 1.434153e+00 1.656316e+00 1.797068e+00 1.978445e+00 3.277463e+00 37.21206 1.0369071
S_0 8.916145e+06 66.4756579 606.5895810 8.914838e+06 8.915746e+06 8.916239e+06 8.916569e+06 8.917209e+06 83.26541 1.0080175
E_0 1.686677e+03 15.9085126 326.3754450 8.965301e+02 1.486263e+03 1.766709e+03 1.949267e+03 2.084575e+03 420.89697 1.0034477
I_{S_0} 2.200000e+01 NA 0.0000000 2.200000e+01 2.200000e+01 2.200000e+01 2.200000e+01 2.200000e+01 NA NA
I_{A_0} 7.993784e+02 72.3470929 554.4589875 1.056505e+02 3.532306e+02 6.434443e+02 1.145775e+03 2.035521e+03 58.73496 1.0048948
H_0 0.000000e+00 NA 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 NA NA
R_0 0.000000e+00 NA 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 NA NA
D_0 0.000000e+00 NA 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 NA NA
ICS_0 7.400000e+01 NA 0.0000000 7.400000e+01 7.400000e+01 7.400000e+01 7.400000e+01 7.400000e+01 NA NA

Lockdown-Vaccination-Treatment-Model

\[\begin{align} S'(t) &= \mu \bar{N} -\frac{\beta_S I_S+\beta_AI_A}{\bar{N}}S -(\mu+\lambda_V)S +\delta_V V \\ E'(t) &= \frac{\beta_S I_S + \beta_A I_A}{\bar{N}} S + \epsilon \frac{\beta_S I_S + \beta_A I_A}{\bar{N}}V - (\mu+\delta_E) E \nonumber \\ I'_S(t) & = p \delta_E E -(\mu + \mu_S + \alpha_S + \lambda_T) I_S +(1 - q) \alpha_T T \\ I'_A(t) & = (1-p) \delta_E E - (\mu + \mu_A + \alpha_A) I_A \\ R'(t) & = \alpha_S I_S + \alpha_A I_A + q \alpha_T T - \mu R \\ D'(t)& = \mu_S I_S + \mu_A I_A \\ V'(t)& = \lambda_V S - \epsilon \frac{\beta_S I_S + \beta_A I_A}{\bar{N}} V - (\mu + \delta_V) V \\ T'(t) &= \lambda_T I_S -(\mu + \alpha_T) T \\ \bar{N}(t) &= S(t) + E(t) + I_S(t) + I_A(t) + R(t) + V(t) + T(t) \end{align}\]

Model from March 24 July 14

state Stan Initial Condition Stan Postulated
\(S\) y[1] \(S_0\) y_init[1] \(N - (I_{S_0} + I_{A_0} + E_0)\)
\(E\) y[2] \(E_0\) y_init[2] Unif(0, 10)
\(I_s\) y[3] \(I_{s_0}\) y_init[3] Unif(0, 3)
\(I_a\) y[4] \(I_{a_0}\) y_init[4] Unif(0, 10)
\(R\) y[5] \(R_0\) y_init[5] Estimated
\(T\) y[5] \(R_0\) y_init[5] 0
\(V\) y[5] \(R_0\) y_init[5] 0
\(I_s^{[cumulative]}\) y[6] \(I_{s_0}^{[cumulative]}\) y_init[6] \(I_{s_0}\)

March 24 –Jun 01 parameters estimation

stan notation symbol Prior Fixed
theta[1] \(\beta_A\) Normal(1, 0.3)
theta[2] \(\beta_S\) Normal(1, 0.3)
theta[3] \(p\) Unif(.3, 0.8)
theta[3] \(q\) Unif(0.4, 0.7)
theta[4] \(\epsilon\) 0.0
theta[5] \(\alpha_S\) Gamma(10, 100)
theta[6] \(\alpha_A\) Gamma(10, 50)
theta[7] \(\alpha_T\) 0.0
theta[8] \(\lambda_V\) 0.0
theta[9] \(\lambda_T\) 0.0
theta[10] \(\mu_S\)
——— \(\mu_A\) 0.0
——— \(\mu\) ——– 3.913894e-05
theta[11] \(\delta_V\) 0.0
theta[12] \(\delta_E\) Gamma(10, 50) 0.1960784
——— \(E_0\) Unif(0, 10)
——— \(I_A\) Unif(0, 10)
——— \(I_S\) Unif(0, 10)
——— \(I_A\) Unif(0, 10)

Stan ODE notation

variable Stan Initial Condition Stan
\(S\) y[1] \(S_0\) y_init[1] \(N - (I{s_0} + I_{a_0} + E_0)\)
\(E\) y[2] \(E_0\) y_init[2] Unif(0, 10)
\(I_s\) y[3] \(I_{s_0}\) y_init[3] Unif(0, 3)
\(I_a\) y[4] \(I_{a_0}\) y_init[4] Unif(0, 10)
\(R\) y[5] \(R_0\) y_init[5] 0
\(D\) y[6] \(D_0\) y_init[6] 0
\(V\) y[7] \(V_0\) y_init[7] 0
\(T\) y[8] \(V_0\) y_init[8] 0
\(I_s^{[cumulative]}\) y[9] \(I_{s_0}^{[cumulative]}\) y_init[9] \(I_{s_0}\)

Stan MODEL

## functions {
##     real[] seirvt(real t,  // time
##              real[] y,
##              real[] theta,
##              real[] x_r,
##              int[] x_i) {
## //
##     real dy_dt[8];
##     real s = y[1];
##     real e = y[2];
##     real i_s = y[3];
##     real i_a = y[4];
##     real h = y[5];
##     real r = y[6];
##     real d = y[7];
##     real n_star;
## //
##     real beta_s = theta[1];
##     real beta_a = theta[2];
##     real kappa = theta[3];
##     real p = theta[4];
##     real delta_h = theta[5];
##     real mu_i_s = theta[6];
##     real mu_h = theta[7];
##     real gamma_s = theta[8];
##     real gamma_a = theta[9];
##     real gamma_h = theta[10];
## //
##     real force_infection;
##     real mu = 3.913894e-05;
##     n_star = s + e + i_s + i_a + h + r;
## 
##     force_infection = (beta_s * i_s + beta_a * i_a) / n_star;
##     dy_dt[1] = mu * n_star - force_infection * s - mu * s;
##     dy_dt[2] = force_infection * s - (mu + kappa) * e;
##     dy_dt[3] = p * kappa * e - (gamma_s + mu_i_s + delta_h) * i_s;
##     dy_dt[4] = (1.0 - p) * kappa * e - (gamma_a + mu) * i_a;
##     dy_dt[5] = delta_h * i_s  - (gamma_h + mu_h + mu) * h;
##     dy_dt[6] = gamma_s * i_s + gamma_a * i_a  + gamma_h * h - mu * r;
##     dy_dt[7] = mu_i_s * i_s + mu_h * h;
##     dy_dt[8] = p * kappa * i_s;
##     return dy_dt;
##   }
## }
## data {
##     int<lower = 1> n_obs;       // number of days observed
##     int<lower = 1> n_theta;     // number of model parameters
##     int<lower = 1> n_difeq;     // number of differential equations
##     int<lower = 1> n_pop;       // population
##     int y[n_obs];               // data, total number of infected
##     real t0;                    // initial time point (zero)
##     real ts[n_obs];             // time points observed
## }
## 
## transformed data {
##   real x_r[0];
##   int x_i[0];
## }
## parameters {
##     real <lower = 5e-6> theta[n_theta];
##     // real <lower = 0, upper = 1> E0;
##     // real <lower = 0, upper=1> IA0;
##     real <lower = 0, upper = n_pop> E0;
##     real <lower = 0, upper = n_pop> IA0;
## }
## transformed parameters{
##     real y_hat[n_obs, n_difeq]; // solution from the ODE solver
##     real y_init[n_difeq];
##     // initial conditions
##     // real IS0 = 8.297217e-06;
##     real IS0 = 22.0;
##     real H0 = 0.0;
##     real R0 = 0.0;
##     real D0 = 0.0;
##     real CIS0 = 74.0;
##     y_init[1] = n_pop - (IS0 + IA0 + E0);
##     y_init[2] = E0;
##     y_init[3] = IS0;
##     y_init[4] = IA0;
##     y_init[5] = H0;
##     y_init[6] = R0;
##     y_init[7] = D0;
##     y_init[8] = CIS0;
##     y_hat = integrate_ode_rk45(seirvt, y_init, t0, ts, theta, x_r, x_i);
## }
## model {
##     real lambda[n_obs];             // poisson parameter
##                                     // priors
##     theta[1] ~ normal(1.0, 0.1);    // beta_s
##     theta[2] ~ normal(1.0, 0.1);    // beta_a
##     theta[3] ~ gamma(10, 40);       // kappa
##     theta[4] ~ uniform(0.1, 0.5);   // p
##     theta[5] ~ gamma(10, 40);       // delta_h
##     theta[6] ~ exponential(33.33333);   // mu_i_s
##     theta[7] ~ exponential(25);   // mu_h
##     theta[8] ~ gamma(10, 100);      // gamma_s
##     theta[9] ~ gamma(10, 50);       // gamma_a
##     theta[10] ~ gamma(10, 250);      // gamma_h
## //
##     // E0 ~ uniform(1.121246e-07, 3.363737e-05);
##     // IA0 ~ uniform(7.848719e-06, 3.363737e-05);
##     E0 ~ uniform(74, 2100);
##     IA0 ~ uniform(74, 2100);
## //    likelihood
##     for (i in 1:n_obs){
##        //lambda[i] = y_hat[i, 8] * n_pop;
##         lambda[i] = y_hat[i, 8];
##     }
##     y ~ poisson(lambda);
## }
## generated quantities {
##     real R_0;      // Basic reproduction number
##     real mu = 3.913894e-05;
##     real beta_s = theta[1];
##     real beta_a = theta[2];
##     real p = theta[4];
##     real gamma_s = theta[4];
##     real gamma_a = theta[5];
##     real kappa = theta[6];
##     R_0 = p * beta_s * kappa / ((gamma_s + mu) * (kappa + mu))
##         + (1.0 - p) * beta_a * kappa / ((gamma_a + mu) * (kappa + mu));
## }

Observation Model

\[\begin{align} Y_t & \sim \text{Poisson}(\lambda_t), \qquad \\ \lambda_t & = \int_0^t p \kappa E \\ & p \sim \text{Uniform(0.3 0.8)} \\ & \kappa \sim \text{Gamma(10, 52)} \end{align}\]